Draft version March 6, 2009 

Preprint typeset using WT^^ style cmulatcapj v. 10/09/06 



o 
o 



O 

u 

6 



> 

m 



m 
o 

o 



X 



STOPPING COOLING FLOWS WITH COSMIC RAY FEEDBACK 

William G. Mathews^ 

Draft version March 6, 2009 

ABSTRACT 

Multi-Gyr two-dimensional calculations describe the gasdynamical evolution of hot gas in the Virgo 
cluster resulting from intermittent cavities formed with cosmic rays. Without cosmic rays, the gas 
evolves into a cooling flow, depositing about 85 solar masses per year of cold gas in the cluster core 
- such uninhibited cooling conflicts with X-ray spectra and many other observations. When cosmic 
rays are produced or deposited 10 kpc from the cluster center in bursts of about 10^^ ergs lasting 20 
Myrs and spaced at intervals of 200 Myrs, the central cooling rate is greatly reduced to M w 0.1 - 
1 solar masses per year, consistent with observations. After cosmic rays diffuse through the cavity 
walls, the ambient gas density is reduced and is buoyantly transported 30-70 kpc out into the cluster. 
Cosmic rays do not directly heat the gas and the modest shock heating around young cavities is offset 
by global cooling as the cluster gas expands. After several Gyrs the hot gas density and temperature 
profiles remain similar to those observed, provided the time-averaged cosmic ray luminosity is about 
Lcr = 2.7 X lO"*"^ erg s^^, approximately equal to the bolometric cooling rate Lx within only ~ 56kpc. 
If an appreciable fraction of the relativistic cosmic rays are protons, gamma rays produced by pion 
decay following inelastic p-p collisions may be detected with the Fermi Gamma Ray Telescope. 

Subject headings: X-rays: galaxies: clusters — galaxies: clusters: general — cooling flows — cosmic 
rays — gamma rays: theory 



1. INTRODUCTION 

X-ray spectra of hot virialized gas in galaxy groups 
and clusters have firmly established that the rate that 
gas cools to low temperatures is much lower than that 
expected from uninhibited radiative cooling in subsoni- 
cally inflowing cooling flows (e.g. Peterson et al. 2001). 
As a result of this realization, many theoretical and com- 
putational studies have sought to prevent cooling inflows 
by introducing a gas heating mechanism usually exploit- 
ing the phenomenal accretion energy released by massive 
black holes in the cores of group and cluster-centered el- 
liptical galaxies (e.g. McNamara & Nulsen 2007). In 
spite of this almost unlimited source of accretion energy, 
no specific mechanism or mode of gas heating has been 
accepted that is successful for long cosmic times and for 
cooling flows of all scales: galactic, group and cluster. 

In their recent review of the cooling flow puzzle McNa- 
mara & Nulsen (2007) (MN07) emphasize two types of 
AGN-related heating mechanisms: PdV work done when 
cavities in the hot gas are inflated by cosmic rays (CRs) 
and energy-dissipating shock (or other) waves that de- 
liver energy from the central energy source throughout 
the cluster gas. CRs, often observed inside cluster cavi- 
ties from radio synchrotron emission, are thought to be 
deposited or produced by non-thermal jets that project 
out from the central AGN. However, our recent study of 
the energetics of cluster cavities inflated by CRs (Math- 
ews & Brighenti 2008b = MB08) shows that the (shock) 
heating produced as young cavities expand due to locally 
enhanced CRs is more than compensated by a global 
cooling as the entire cluster gas expands more or less 
permanently to accommodate the partial pressure of the 
CRs as they diffuse into the cluster gas. Most of the 
PdV work expended as cavities form is stored in poten- 

^ UCO/Lick Observatory, Dept. of Astronomy and Astro- 
physics, University of California, Santa Cruz, CA 95064 



tial energy of displaced gas that moves out in the clus- 
ter potential. The dissipative heating beneath buoyant 
cavities discussed in MN07 appears to be absent in the 
cavities computed in MB08, and this is consistent with 
our earlier estimate (Mathews et al. 2006). The total 
(bulk and turbulent) kinetic energy is always very small 
and cannot be an important source of dissipative heating 
(MB08). 

The second currently favored heating mechanism, large 
scale wave dissipation, is also problematic. Since the hot 
gas density profiles in groups and clusters is generally 
flatter than p cx r^^, successive shocks preferentially heat 
gas near the cluster center until it exceeds observed tem- 
peratures (Mathews, Faltenbacher & Brighenti 2006). 
Over-heating the central gas is a persistent difficulty 
with previous theoretical attempts to resolve the cool- 
ing flow problem by using central AGN heating scenar- 
ios (Brighenti & Mathews 2002b, 2003) since the cluster 
gas closest to the source of AGN energy typically has 
the lowest temperatures observed in the hot cluster gas. 
Finely-tuned heating scenarios in which the gas inflow is 
assumed to be essentially stopped are inconsistent with 
secular changes in the gas that evolve without limit. For 
example metal enrichment by Type la supernovae in the 
central galaxy would eventually increase the local hot gas 
metal abundances far above those observed. 

As an alternative solution to the cooling flow problem, 
we proposed that gas receiving AGN energy is buoyantly 
transported outward. In these mass-circulating flows gas 
flows in both radial directions simultaneously. The com- 
bined flows maintain a cooling flow appearance similar 
to observed X-ray images (Mathews et al. 2003, 2004; 
Brighenti & Mathews 2006). 

In this Letter we describe the evolution of virialized 
hot gas in the Virgo cluster in which successive cavities 
are formed with CRs that buoyantly transfer gas from 
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near the cluster core to large radii where its radiative 
cooling time is very much longer. The average net gas 
mass that flows past every cluster radius is very small 
although significant inflows and outflows occur locally in 
the gas. In buoyant outflows generated by CRs, the clus- 
ter gas need not be directly heated by the CRs, although 
some modest heating results from the dissipation of weak 
shocks that propagate away from newly formed cavities 
(Mathews & Brighenti 2007; MB08). UnUke most AGN 
heating scenarios, buoyant mass-circulating flows gener- 
ated by CRs can preserve the observed cluster tempera- 
ture and density profiles for at least several Gyrs. This 
resolution of the cooling flow problem, using mass trans- 
fer by CR buoyancy rather than heating, is attractive 
because almost all group and cluster-centered galaxies 
are observed to produce CRs in jets, radio lobes, or in 
central non-thermal radio sources. 

2. EQUATIONS AND PROCEDURE 

The procedure and assumptions employed here are 
nearly identical to those described in detail by MB08 
to which the reader is referred. The flow equations are 

|+V.pu = (1) 
pf^ + (u-V)u') =-V(P + Pc)-pg (2) 
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where artiflcial viscosity terms are suppressed. Gas and 
CR pressures are related to their corresponding energy 
densities by P = (7 — l)e and Pc = (7c — l)eo respectively 

where for simplicity we assume 7 = 5/3 and 7^ = 4/3 for 
fully relativistic CRs. The energy spectrum of the CR 
particles N(E) is not considered in detail so the CRs 
are characterized only by their integrated energy density 
Cc oc / EN{E)dE. Furthermore, the CR particles are 
unspecified and can consist of electrons or protons in any 
combination. 

Unlike the discussion in MBG8, we now include ra- 
diative cooling. We adopt the radiative cooling coef- 
ficient nionsne^sdiT, z) erg cm~^ s~^ from Sutherland 
and Dopita (1993) which, as expressed in equation 3, 
{p/mpY\ erg cm^'^ s~^, requires that A — 1.1[(4 — 
2>ii){2 + ^)/25/x^]Asd where /x = 0.61 is the molecular 
weight and the coefficient is Uions/np = 1.1. We do 
not compute the metallicity in the gas throughout the 
fiow but simply adopt a uniform metallicity near so- 
lar z = O.75z0 since most of the gas that cools has 
been metal- enriched in the cluster core. The computed 
flows are not sensitive to the gas metallicity. Computa- 
tional timesteps are required to be less than the constant- 
pressure radiative cooling time in all zones. Since some 
cooling occurs in off-center zones near the z-axis, it is use- 
ful to remove cooled gas from the computational grid. If 
a grid zone cools below Tmin = 5 x 10^ K, gas is removed 
at constant pressure by reducing the zone density until 
the temperature resets to Tmax = 5 x 10^ K. Gas removal 
at constant pressure from a zone preserves the pressure 
gradients relative to nearby zones so the gas flow to or 



from neighboring zones is unaffected. Moreover, the to- 
tal thermal energy in the zone is also preserved so thor- 
oughly cooled gas removed in this way does not alter the 
total thermal energy of uncooled gas. We assume that 
CRs are not removed with the cooled gas. 

Gas is accelerated by pressure gradients in both the 
thermal gas and CRs. The two fluids are assumed to 
be coupled by the presence of small magnetic fields in 
the cluster gas. This coupling can be efficient even when 
the energy density in the magnetic field is much less than 
that in the thermal gas, which is consistent with the small 
fields observed in clusters, 1 — 10/iG (Govoni & Fcrctti 
2004). Consequently, magnetic pressure and stresses are 
not explicitly considered. 

Equation 4 above describes the advection of CRs by 
locally flowing gas and the diffusion of CRs through this 
gas. Little is known about the diffusion coefficient k, 
but it is likely that it varies inversely with the gas den- 
sity since the magnetic field is probably larger in denser 
gas. As in MB08 we consider a range of CR diffusion 
coefficients that depend on a a single density parameter 
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Here we consider rieo = 6 x 10"'^ and 6 x 10~^ cm^^. The 
cluster gas dynamics are not very sensitive to rieo or to 
other K(r7,e) variations that have been explored (MB08). 

The flow equations are solved in 2D cylindrical coor- 
dinates [R, z) using a ZEUS-likc code (Stone & Norman 
1992) but with the additional CR equation. As in MB08 
the computational grid consists of 100 equally-spaced 
zones out to 50 kpc plus an additional 100 zones in both 
coordinates that increase logarithmically out to ~ 1 Mpc, 
filling a large hemisphere. Also following MB08, a spher- 
ical gravitational potential for the Virgo cluster is found 
by assuming that the gas is in hydrostatic equilibrium 
with the temperature and modified density profiles Tc{r) 
and PrnGir) observed by Ghizzardi et al. (2004). The 
density profile was found by integrating the equation of 
hydrostatic equilibrium using Tcir) and the total grav- 
itational acceleration from a 3.2 x 10^ Mq black hole, 
the M87 stellar mass p*{r) from Brighenti & Mathews 
(2002a), and an NFW profile with M^ir = 1.30 x 10^^ 
Adr.) and concentration 8.60. The resulting modified ob- 
served density profile PmG{T) agrees very well with that 
of Ghizzardi et al. pa {r) for r > 1 kpc but avoids an un- 
physical maximum in g{r) near ~ 4 kpc if g{r) is found 
from Tcir) and pai^)- All flow calculations begin with 
Tcir) and Pmcir)- 

CRs are introduced at intervals of tcyc — 2 x 10* yrs 



in bursts lasting tc 



2 X 10' yrs. These timescales 



are roughly consistent with the lifetimes of visible X-ray 
cavities and radio synchrotron electrons. If the CR lu- 
minosity during injection times is Ec, the time-averaged 
CR luminosity is 



{-^cr} — 2(tc(iv /^^cycle) -^c- 



(5) 



assuming that CRs are injected into both cluster hemi- 
spheres by symmetric double jets. For comparison with 
the cavity calculations in MB08, all cavities are assumed 
to form in a source region 10 kpc along the z— axis de- 
scribed by a normalized Gaussian profile. The CR source 
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term in equation 4 is therefore 

e-((^-^cav)/rsf 
S^ = E^T{t) -— ergcm-3s-i (6) 

where r is the distance from the origin, r, = 2 kpc and 
T{t) is unity during injection times and zero otherwise. 
All calculations described here continue for 3 Gyrs. For 
each assumed K(neo) we vary E^, seeking that value for 
which the gas temperature and density profiles continue 
to agree best with profiles observed in Virgo. Except 
for radiative cooling and a variable Ec, the parameters 
chosen here are identical to those discussed for single 
cavities by MB08. 

3. RESULTS 

The left column of Figure 1 illustrates the evolution 
of Virgo cluster gas as a pure radiatively cooling flow 
without CRs. The long-dashed lines show the radial gas 
density, temperature and pressure profiles after evolv- 
ing for 3 Gyrs from the observed profiles. The evolved 
profiles are determined by ordering p, T and P in each 
computational zone by increasing distance from the ori- 
gin r = {R^ + z^)^/^, then averaging in appropriate 
bins Ar. After 3 Gyrs, p{r) and T{r) have evolved 
away from the initial (observed) profiles (solid lines). 
The ratio of gas entropy to its initial value (subscript 
0), log(5/So) = log(T/To) - (2/3)log(/9/po), decreases 
monotonically toward the cluster center where the gas 
cools in the central 1.5 kpc. The discrepancy in inte- 
grated mass of (uncooled) gas between the initial and 
evolved density profiles in the uppermost panel, AM = 
M(r,3 Gyr) - M(r,0), peaks at ~ 15 kpc^. Between 
2 and 3 Gyrs the quasi-steady cooling rate in the core 
is (M) = 86 Mq yr~^ for the entire Virgo cluster (both 
hemispheres)^. 

Panels in the central column of Figure 1 show the corre- 
sponding evolution with CRs having a low CR diffusion 
coefficient n with rieo = 6 x 10~^ cm~^. In this solu- 
tion after 3 Gyrs the gas density and temperature re- 
main close to the initial, observed profiles. The entropy 
in 25 r 100 kpc is lowered by the accumulated buoy- 
ant arrival of gas from the cluster core (cf. MB08) that 
lost entropy by radiation. The modest entropy maximum 
near 10 kpc is probably due to heating by successive weak 
cavity shocks at this radius. The total mass of hot gas 
in each spherical annulus changes very little during the 
evolution, as expected from the small central mass cool- 
ing rate between 2 and 3 Gyrs, {M) = 0.13 Mq yr^^ (for 
both hemispheres). The cooling inflow has been almost 
completely balanced by the outward circulation of gas 
from the cluster core to large radii driven by CR buoy- 
ancy. The bottom panel shows that the CR pressure Pc 
after 3 Gyrs is about an order of magnitude less than the 
gas pressure within 100 kpc, but Pc/ P decreases sharply 
beyond 100 kpc. In this solution the total (spherical) 
time-averaged CR luminosity is (Lcr) = 2.70 x 10''^ erg 
s~^. The total CR energy delivered to each cavity is 
about E^av = 8.2 x 10^^ ergs. 

^ This suggests the approximate site to initiate buoyant outflow. 

* The mass coohng rate, found from tfie total mass of eooled 
gas removed during this time interval, is insensitive to T^in arid 
Tmax- For example, (M) = 85 Mq yr^^ with Tmin = 3 X 10^ K 
and Tmax = 10« K. 



The right column of panels in Figure 1 shows the iden- 
tical evolution but with a much larger CR diffusion coef- 
ficient (rieo = 6 X 10~^ cm~^). As before, the gas density 
and temperature profiles remain essentially unchanged 
after 3 Gyrs with a total time-averaged CR luminosity 
(Lcr) = 2.60 X 10''^ erg s~^. The total gas cooling rate 
between 2 and 3 Gyrs, (M) = 1.15 Mq yr~^, is only 
~ 0.01 of the cooling flow rate without CRs. 

4. DISCUSSION AND CONCLUSIONS 

CR buoyancy efficiently arrests cooling inflows with- 
out requiring ad hoc heating or infusions of hot, non- 
relativistic gas in the cavities. In the two solutions with 
CRs gas cools at the center and along the z-axis out to 
about 3-5 kpc where the free fall time to the center is 
< 10'' yrs. If all this cooled gas were accreted by the cen- 
tral black hole, the efficiency required for accretion en- 
ergy to produce enough CRs to maintain the observed gas 
profiles is rjcr = {Lcr)/{M)c^. We find jjcr = 3.6 x 10~^ 
and 3.9 x lO"** respectively for the low and high k solu- 
tions. But rjcr is probably underestimated since some of 
the cooled, dusty gas in group-centered elliptical galaxies 
is heated and buoyantly relocated far out into the hot gas 
(Temi, Brighenti, & Mathews 2007). The total CR lumi- 
nosity required to (nearly) stop the Virgo cooling flow, 
(Lcr) = 2.60 X lO"'^ erg s~^, is equal to the (computed) 
bolometric X-ray luminosity within only ^ 56kpc. 

The Fermi Gamma Ray Telescope may confirm that 
cooling inflows in Virgo and other nearby clusters are 
stopped with CR buoyancy. Gamma rays are created by 
neutral pion decays following inelastic p-p interactions 
between CRs and the thermal plasma. The integrated 
gamma ray flux (cm^^ s^^) above energy E^ „iin is 

F^{E^ > E^^rnin) = / dV ( dE^q^{E^,r) 

47ra Jvirgo J 

(7) 

where d = 16 Mpc is the distance to Virgo and the spe- 
cific photon emissivity, qj{Ej,r) cm~^ GeV~^, is 
described by Pfrommer & Ensslin (2004) and Ando & 
Nagai (2008). Assume for simplicity that all CRs are 
relativistic protons having a power law momentum dis- 
tribution f{pp, r) oc with a — 2.3. Then after 3 Gyrs 
the integrated gamma ray photon flux from Virgo above 
E.y,min = 2 Gev isF^ = 2x IQ-^ or 3 x lO'^ cm'^ s'^ 
respectively for the low and high k solutions in Figure 
1. These fluxes are near the detection threshold of the 
Fermi telescope, ~ 10^^ cm^^ s^^ estimated by Ando 
& Nagai (2008) for extended sources. If half of the CRs 
are electrons, evidence of the AGN CRs that we require 
in Virgo may still be marginally visible with the Fermi 
Telescope. Additional CRs in Virgo may have been pro- 
duced in shocks and turbulence accompanying mergers 
or in the accretion shock when baryons flrst fell into the 
Virgo potential (Loeb & Waxman 2000), but these are 
not considered here. 

For both low and high k solutions in Figure 1 the CR 
partial pressure Pc/{P + Pc) increases toward the cluster 
center. This may seem surprising because all CRs are 
injected at 10 kpc and diffusion into the dense, central 
gas is slow because k decreases. This feature can be un- 
derstood from the reaction of X-ray cavities on the clus- 
ter gas. As a young cavity buoyantly moves away from 
the cluster center, relatively cooler nearby gas concen- 
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trates underneath the cavity, then flows rapidly outward 
through the cavity along the z-axis (MB08). These jets 
or filaments of denser thermal gas move out beyond 10 
kpc then later fall back under gravity into the cluster 
core (e.g. Mathews & Brighenti 2008a). It is likely that 
the central concentration of CRs occurs because they are 
advected in with each infalling filament, then become 
trapped by the relatively small k in the dense central 
gas. CRs arc a likely source of non-thermal pressure that 
is occasionally observed in the cores of hot galactic and 
cluster gas. 

Figure 2 shows that the distribution of CRs in in Virgo 
after 3 Gyrs is concentrated in the core and along the 
^;-axis. In a recent study of warm and cold gaseous fil- 
aments in the Perseus Cluster, Ferland et al. (2008a,b) 
find that the observed spectrum can be understood if 
some of the the emission lines are excited by collisions 
with CRs. The enhanced CR distribution in Figure 2 
along the z-axis provides a natural source of CRs just 
where the cool radial filaments form. Although we do 
not compute the passively evolving magnetic field, it is 
expected to be stronger in denser thermal filaments and 
longitudinally oriented along the cool filaments due to 
the radial expansion of the filament fiow along the z-axis. 

The 2D cylindrical nature of our solution requires that 
all cavities and CR activity occur along the z-axis. The 
relatively small number of double-double radio sources 
(Schoenmakers et al. 2000) indicates that successive CR 
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jets can move along the same direction, as assumed here. 
Remarkably, the X-ray thermal filament in Virgo (For- 
man et al. 2007) that formed by a cavity ~ 10^ years 
ago (Mathews & Brighenti 2008a) and its associated ra- 
dio lobe lie in a different direction than the younger non- 
thermal jet in the central galaxy M87. Another indica- 
tion that CR energy flows out in many different direc- 
tions is the azimuthal variation of successive generations 
of X-ray cavities in Perseus and other clusters. 

It is expected that the mass outflow due to CR buoy- 
ancy in Figure 1 also transports iron produced by Type la 
supernova inside the central galaxy into a quasi-spherical 
region 50-100 kpc in radius as observed in many groups 
and clusters (De Grandi et al. 2004). The outflows de- 
scribed here can also carry short-lived dust out to 5-10 
kpc into the hot gas where its far-infrared emission is ob- 
served in many nearby galaxy groups (Temi, Brighenti, 
& Mathews 2007). Radial entropy variations qualita- 
tively similar to those in Figure 1 have been observed in 
galaxy groups (Gastaldello et al. 2007; Sun et al. 2008). 
Finally, the outward mass transport can create an in- 
flection in the density slope as in Figure 1, resembling 
the double-structured character of many density profiles 
observed. 

Studies of the evolution of hot cluster gas at UC Santa 
Cruz are supported by NASA and NSF grants for which 
we are very grateful. 
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Fig. 1. — Radial distribution of gas and CRs in Virgo after 3 Gyrs of intermittent CR injection. Columns from left to right show a 
cooling flow without CRs, the same flow with weakly then strongly diffusing CRs. Top row: initial (observed) gas density {solid lines) and 
density after 3 Gyrs {long dashed lines). Central row: initial (observed) gas temperature {solid lines), and after 3 Gyrs {long dashed lines) 
(in units of lO'^ K) and the change log S/ So of the gas entropy after 3 Gyrs {dotted lines). Bottom row: pressure of gas {long dashed lines) 
and CRs {short dashed lines) after 3 Gyrs (units of 10"^" dynes cm~^). 
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Fig. 2. — Distributions after 3 Gyrs of CR energy density ec(R, z) in units of 10 erg cm ^ (solid lines) and cooled mass in 10* Mq 
{dashed lines). 



